clear all
clc
%Initial parameters

%State dependent utility pi
pi_g=0;
pi_b=-1;

%---------Five exogenuos variables---------
mu=[0.1 0.5, 0.9]; %probability of good state
Num_mu=size(mu,2);

lambda=[0.01,0.02]; %high and low type of Shannon cost
Num_lambda=size(lambda,2);

MGrid=199;
m=(1/MGrid):((pi_g-pi_b)/MGrid):(pi_g-pi_b-1/MGrid); %material incentive

sigma=zeros(Num_lambda,MGrid-1,Num_mu); %optimal sigma
prtc=zeros(Num_lambda,MGrid-1,Num_mu); %ex ante participation probability
decision=zeros(Num_lambda,MGrid-1,Num_mu); %record decisions
gammaA=zeros(Num_lambda,MGrid-1,Num_mu); %ex ante participation probability
gammaR=zeros(Num_lambda,MGrid-1,Num_mu); %ex ante participation probability

   
 

for k=1:1:Num_mu
    for i=1:1:(MGrid-1)
        for j=1:1:Num_lambda
            gamma = -mu(k)/(1-mu(k))*(pi_g+m(i))/(pi_b+m(i));
            a = (log(gamma))^2;
            b = log(-2*pi*lambda(j)^2/mu(k)/(1-mu(k))/(pi_g+m(i))/(pi_b+m(i)));
            c = 1/4;
            delta = b^2 - 4*a*c;
            Uprior = mu(k)*(pi_g+m(i))+(1-mu(k))*(pi_b+m(i));
            if b<0 && delta>0
                x1 = sqrt((-b-sqrt(delta))/(2*a));
                U = mu(k)*normcdf((0.5+x1^2*log(gamma))/x1)*(pi_g+m(i))+(1-mu(k))*(1-normcdf((0.5-x1^2*log(gamma))/x1))*(pi_b+m(i))-lambda(j)/x1;
                p_s = mu(k)*normcdf((0.5+x1^2*log(gamma))/x1)+(1-mu(k))*(1-normcdf((0.5-x1^2*log(gamma))/x1));

                pG = normcdf((0.5+x1^2*log(gamma))/x1);
                pB = (1-normcdf((0.5-x1^2*log(gamma))/x1));
                

                if gamma<=1 && U>=0
                    sigma(j,i,k) = x1;
                    prtc(j,i,k) = p_s;
                    decision(j,i,k) = 1; %1 for choose using sigma_star, ex ante not participate
                    gammaA(j,i,k) = mu(k)*pG ./ (mu(k)*pG + (1 - mu(k))*pB);
                    gammaR(j,i,k) = mu(k)*(1 - pG) ./ (mu(k)*(1 - pG) + (1 - mu(k))*(1 - pB));
                elseif gamma<=1 && U<0
                    sigma(j,i,k) = NaN;
                    prtc(j,i,k) = 0;
                    decision(j,i,k) = 5; %5 for not participate, i.e. sigma_star is infinity
                    gammaA(j,i,k) = NaN;
                    gammaR(j,i,k) = NaN;
                elseif U>=Uprior
                    sigma(j,i,k) = x1;
                    prtc(j,i,k) = p_s;
                    decision(j,i,k) = 2; %2 for choose using sigma_star, ex ante participate
                    gammaA(j,i,k) = mu(k)*pG ./ (mu(k)*pG + (1 - mu(k))*pB);
                    gammaR(j,i,k) = mu(k)*(1 - pG) ./ (mu(k)*(1 - pG) + (1 - mu(k))*(1 - pB));
                else
                    sigma(j,i,k) = NaN;
                    prtc(j,i,k) = 1;
                    decision(j,i,k) = 3; %3 for participate without getting the signal
                    gammaA(j,i,k) = NaN;
                    gammaR(j,i,k) = NaN;
                end
            else
                sigma(j,i,k) = NaN;
                prtc(j,i,k) = (gamma>1);
                gammaA(j,i,k) = NaN;
                gammaR(j,i,k) = NaN;
                %decision 0 means no optimal solution for sigma
            end
        end
    end
    
    
 
    frac_htype=prtc(2,:,k)./(prtc(1,:,k)+prtc(2,:,k));
    L=plot(m,prtc(1,:,k),'r',m,prtc(2,:,k),'b--',m,frac_htype,'k');
    set(L(1),'LineWidth',2);
    set(L(2),'LineWidth',2);
    set(L(3),'LineWidth',4);
    xlabel('Incentive m');
    ylabel('Participation probability');
    plotname1=['\lambda = ' num2str(lambda(1))];
    plotname2=['\lambda = ' num2str(lambda(2))];
    legend(plotname1,plotname2,'P(cost = high | participate)');
    if k >= 2
       legend('Location','southeast');
    else
       legend('Location','northwest');
    end
    prior=mu(k)*100;
    filename=['../graphs/supplyNormalEndo' num2str(prior) '.pdf'];
    saveas(gcf,filename); 
    
    L=plot( squeeze(m), squeeze(gammaA(1,:,k))       ,'r' , ...
            squeeze(m), squeeze(gammaA(2,:,k))       ,'b--', ...
            squeeze(m), squeeze(gammaR(1,:,k))       ,'r' , ...
            squeeze(m), squeeze(gammaR(2,:,k))       ,'b--' ...
            );

        
    set(L(1),'LineWidth',2);
    set(L(2),'LineWidth',2);
    set(L(3),'LineWidth',2);
    set(L(4),'LineWidth',2);

    xlabel('Incentive m');
    ylabel('P(s = G | action)');
    plotname1=['\lambda = ' num2str(lambda(1))];
    plotname3=['\lambda = ' num2str(lambda(2))];
    legend(plotname1,plotname3);
    if mu(k) < 0.55
        legend('Location','west');
    else
        legend('Location','east');
    end
    prior=mu(k)*100;
    axis manual 
    axis([0 1 0 1])
    filename=['../graphs/posteriorsNormalEndo' num2str(prior) '.pdf'];
    saveas(gcf,filename); 
    
end



